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We propose a time-dependent many-body approach to study the short-time dynamics of corre- 
lated electrons in quantum transport through nanoscale systems contacted to metallic leads. This 
approach is based on the time-propagation of the Kadanoff-Baym equations for the nonequilibrium 
many-body Green's function of open and interacting systems out of equilibrium. An important 
feature of the method is that it takes full account of electronic correlations and embedding effects 
in the presence of time-dependent external fields, while at the same time satisfying the charge 
conservation law. The method further extends the Meir-Wingreen formula to the time domain for 
initially correlated states. We study the electron dynamics of a correlated quantum wire attached to 
two-dimensional leads exposed to a sudden switch-on of a bias voltage using conserving many-body 
approximations at Hartree-Fock, second Born and GW level. We obtain detailed results for the 
transient currents, dipole moments, spectral functions, charging times, and the many-body screen- 
ing of the quantum wire as well as for the time-dependent density pattern in the leads, and we 
show how the time-dependence of these observables provides a wealth of information on the level 
structure of the quantum wire out of equilibrium. For moderate interaction strenghts the 2B and 
GW results are in excellent agreement at all times. We find that many-body effects beyond the 
Hartree-Fock approximation have a large effect on the qualitative behavior of the system and lead 
to a bias dependent gap closing and quasiparticle broadening, shortening of the transient times and 
washing out of the step features in the current-voltage curves. 

PACS numbers: 72.10.Bg,71.10.-w,73.63.-b,85.30.Mn 



I. INTRODUCTION 

The description of electron transport through 
nanoscale systems contacted to metallic leads is cur- 
rently under intensive investigation especially due to 
the possibility of miniaturizing integrated devices in 
electrical circuitsji Several theoretical methods have 
been proposed to address the steady state properties of 
these systems. 

Ab initio formulations based on Time-Dependent (TD) 
Density Functional Th.eoTy^'^i^^^ (DFT) and Current 
Density Functional Theory^i^ i^°'^^i^^ provide a virtual ex- 
act framework to account for correlation effects both in 
the leads and the device but lack of a systematic route to 
improve the level of the approximations. Ad hoc approx- 
imations have been successfully implemented to describe 
qualitative features of the I/V characteristic of molecular 
junctions in the Coulomb blockade regime J^'i^ii^'i^ More 
sophisticated approximations are, however, needed for, 
e.g., non-resonant tunneling transport through weakly 
coupled molecules 

The possibility of including relevant physical processes 
through an insightful selection of Feynman diagrams is 
the main advantage of Many-Body Perturbation The- 
ory (MBPT) over one-particle schemes. Even though 
computationally more expensive MBPT offers an invalu- 
able tool to quantify the effects of electron correlations 
by analyzing, e.g., the quasi-particle spectra, life-times. 



screened interactions, etc. One of the most remarkable 
advances in the MBPT formulation of electron transport 
was given by Meir and Wingreen who provided an equa- 
tion for the steady state current through a correlated de- 
vice regioEk^ii^ thus generalizing the Landauer formulai^ 
The Meir-Wingreen formula is cast in terms of the in- 
teracting Green's function and self-energy in the device 
region and can be approximated using standard diagram- 
matic techniques. Exploiting Wick's theorem^! a general 
diagram for the self-energy can be written in terms of 
bare Green's functions and interaction lines. Any ap- 
proximation to the self-energy which contains a finite 
number of such diagrams does, however, violate many 
conservation laws. Conserving approximations2^i^i21i^ 
require the resummation of an infinite number of dia- 
grams and are of paramount importance in nonequilib- 
rium problems as they guarantee satisfaction of funda- 
mental conservation laws such as charge conservation. 
Examples of conserving approximations are the Hartree- 
Fock (HF), second Born (2B), GW, T-matrix, and fluctu- 
ation exchange (FLEX) approximations The success 
of the GW approximation^ii^^ in describing spectral fea- 
tures of atoms and molccules^i^ii^^ as well as of interact- 
ing model clusters^^ prompted efforts to implement the 
Meir-Wingreen formula at the GW level in simple molec- 
ular junctions and tight-binding modelsi ^^'^^i'^^'^°'^^i"^^ 

The advantage of using molecular devices in future na- 
noelectronics is, however, not only the miniaturization of 
integrated circuits. Nanodcvices can work at the THcrtz 
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regime and hence perform operations in few picoseconds 
or even faster. Space and time can both be consider- 
ably reduced. Nevertheless, at the sub-picosecond time 
scale stationary steady-state approaches are inadequate 
to extract crucial quantities like, e.g., the switching- or 
charging-time of a molecular diode, and consequently to 
understand how to optimize the device performance. De- 
spite the importance that an increase in the operational 
speed may have in practical applications, the ultrafast 
dynamical response of nanoscale devices is still largely 
unexplored. This paper wants to make a further step 
towards the theoretical modeling of correlated TD quan- 
tum transport. 

Recently several practical schemes have been proposed 
to tackle TD quantum transport problems of noninter- 
acting electrons i^i^i^i^i^ In some of these schemes 
the electron-electron interaction can be included within 
a TDDFT framework^^ and few calculations on the 
transient electron dynamics of molecular junctions have 
been performed at the level of the adiabatic local density 
approximation . ^^1^^'^° Alternatively, approaches based 
on Bohm trajectories^i*^ or on the density matrix renor- 
malization groupS have been put forward to calculate 
TD currents and densities through interacting quantum 
systems. So far, however, no one has extended the dia- 
grammatic MBPT formulation of Meir and Wingrccn to 
the time-domain. As in the steady-state case the MBPT 
formulation allows for including relevant scattering mech- 
anisms via a proper selection of physically meaningful 
Feynman diagrams. The appealing nature of diagram- 
matic expansions renders MBPT an attractive alterna- 
tive to investigate out-of-equilibrium systems. 

In a recent Letter— we proposed a time-dependent 
MBPT formulation of quantum transport which is based 
on the real-time propagation of the Kadanoff-Baym (KB) 
equation o^^'^^i^'^'^^i^^'^^i^^ for open and interacting sys- 
tems. The KB equations are equations of motion for the 
nonequilibrium Green's function from which basic prop- 
erties of the system can be calculated. It is the purpose 
of this paper to give a detailed account of the theoret- 
ical derivation and to extend the numerical analysis to 
quantum wires connected to two-dimensional leads. For 
practical calculations we have implemented the fully self- 
consistent HF, 2B and GW conserving approximations. 
Our results reduce to those of steady-state MBPT imple- 
mentations in the long time limit. Having full access to 
the transient dynamics we are able, however, to extract 
novel information like the switching- and charging-times, 
the time-dependent renormalization of the electronic lev- 
els, the role of initial correlations, the time-dependent 
dipole moments etc. Furthermore, the non-locality in 
time of the 2B and GW self-energies allows us to highlight 
non-trivial memory effects occuring before the steady- 
state is reached. We also wish to emphasize that our 
approach is not limited to DC biases. Arbitrary driving 
fields like AC biases, voltage pulses, pumping fields, etc. 
can be dealt with at the same computational cost. 

The paper is organized as follows. All derivations and 



formulas are given in Section |TT1 We present the class of 
many-body systems that can be studied within our KB 
formulation in Section III Al and derive the equations of 
motion for the nonequilibrium Green's function in the 
device region in Section [ll Bl (see also Appendix [A]) . The 
equations of motion are then used to prove the continuity 
equation for all conserving approximations. Section [II C[ 
and to extend the Meir-Wingrcen formula to the time do- 
main for initially correlated systems. Section [II Dl Using 
an inbedding technique in Section jTlE] we derive the main 
equations to calculate the time-dependent density in the 
leads. In Section IIIII we present the results of our TD 
simulations for a one-dimensional wire connected to two- 
dimensional leads. The Keldysh Green's function, which 
is the basic quantity of the KB approach, of the open 
wire is studied in Section IIII Al showing different time- 
dependent regimes relevant to the subsequent analysis. 
In Section IIII Bl and IIII CI we calculate the TD current 
and dipole moment respectively. We find that the 2B and 
GW results are in excellent agreement at all times and 
can differ substantially from the HF results. We also per- 
form the Fourier analysis of the transient oscillations and 
reveal the underlying out-of-equilibrium electronic struc- 
ture of the open wirci^ The dynamically screened inter- 
action of the GW approximation is investigated in Sec- 
tion |TTTD] with emphasis on the time-scales of retardation 
effects. Section llll El is devoted to the study of the TD re- 
arrangement of the density in the two-dimensional leads 
after the switch-on of an external bias. Such an analy- 
sis permits us to test the validity of a commonly used 
assumption in quantum transport, i.e., that the leads re- 
main in thermal equilibrium. Finally, in Section IIVI we 
draw our main conclusions and future perspectives. 

II. THEORY 

A. The model Hamiltonian 

We consider a class of quantum correlated open sys- 
tems (which we call central regions) coupled to noninter- 
acting reservoirs (which we call leads), see Fig. [TJ The 
Hamiltonian has the general form 

H{t)=Hc{t)+J2Hc.{t)+HT-t^N, (1) 

a 

where Hq, Ha, Ht are the central region, the lead a 
and the tunneling Hamiltonians respectively and N is the 
particle number operator coupled to chemical potential ^. 
We assume that there is no direct coupling between the 
leads. The correlated central region has a Hamiltonian 
of the form 

Hc{i) = '^^dt)dldja + ^ J2 v^3kldld]a'dka'dla, (2) 

ij,(T ijkl 

where i,j label a complete set of one-particle states in the 
central region, a, a' are spin- indices and d) , d are the ere- 
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FIG. 1; Sketch of the transport setup. The correlated cen- 
tral region (C) is coupled to semi-infinite left (L) and right 
(R) tight-binding leads via tunneling Hamiltonians H^q and 
Hca, a = L, R. 



ation and annihilation operators respectively. The one- 
body part of the Hamiltonian hij(t) may have an arbi- 
trary time-dependence, describing, e.g., a gate voltage or 
pumping fields. The two-body part accounts for interac- 
tions between the electrons where Vijki are, for example 
in the case of a molecule, the standard two-electron inte- 
grals of the Coulomb interaction. The lead Hamiltonians 
have the form 



(3) 



where the creation and annihilation operators for the 
leads are denoted by c'^ and c. Here Na — J2i a ^laa^ic^a 
is the operator describing the number of particles in lead 
a. The one-body part of the Hamiltonian hfj describes 
metallic leads and can be calculated using a tight-binding 
representation, or a real-space grid or any other conve- 
nient basis set. We are interested in exposing the leads 
to an external electric field which varies on a time-scale 
much longer than the typical plasmon time-scale. Then, 
the coarse-grained time evolution can be performed as- 
suming a perfect instantaneous screening in the leads and 
the homogeneous time-dependent field Ua{t) can be in- 
terpreted as the sum of the external and the screening 
field, i.e., the applied bias. This effectively means that 
the leads are treated at a Hartree mean field level. We 
finally consider the tunneling Hamiltonian Ht 



(4) 



which describes the coupling of the leads to the interact- 
ing central region. This completes the full description 
of the Hamiltonian of the system. In the next section 
we study the equations of motion for the corresponding 
Green's function. 



Equation of motion for the Keldysh Green's 
function 



We assume the system to be contacted and in equi- 
librium at inverse temperature (3 before time t = to and 
described by Hamiltonian Hq. For times t > to the sys- 
tem is driven out of equilibrium by an external bias and 
we aim to study the time-evolution of the electron den- 
sity, current, etc.. In order to describe the electron dy- 
namics in this system we use Keldysh Green's function 
theory (for a review see Reflsol) which allows us to include 
many-body effects in a diagrammatic way. The Keldysh 
Green's function is defined as the expectation value of 
the contour-ordered product 



Grs{z, z) 



Tr 



Tr 



(5) 



where a and are either lead or central region operators 
and the indices r and s are collective indices for position 
and spin. The variable z is a time contour variable that 
specifies the location of the operators on the time con- 
tour. The operator T orders the operators along the 
Keldysh contour displayed in Fig. [21 consisting of two 
real time branches and the imaginary track running from 
to to to — ip. In the definition of the Green's function 
the trace is taken with respect to the many-body states 
of the system. 

All time-dependent one-particle properties can be calcu- 
lated from Q. For instance, the time-dependent density 
matrix is given as 



nrs{t) = -i0rs{t-,t+). 



(6) 



where the times t± lie on the lower/upper branch of the 
contour. The equations of motion for the Green's func- 
tion of the full system can be easily derived from the 
definition Eq. ^ and read 

dz-E^^iz,z)giz,z'), (7) 



(8) 



-idz'G{z,z') = S{z,z')l + g{z,z')H{z') 

+ J dzg{z,z)i:^''{z,z), 



where S is the many-body self-energy, H(z) is the 
matrix representation of the one-body part of the full 
Hamiltonian and the integration is performed over the 
Keldysh-contour. This equation of motion needs to be 
solved with the boundary conditions^i^ 



g{h,z') 
g{z,to) 



-g{to - iP, z ), 
-g{z,to - if3), 



(9) 



which follow directly from the definition of the Green's 
function Eq. Explicitly, the one-body Hamiltonian 
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FIG. 2; Keldysh contour 7. Times on the upper/lower branch 
are specified with the subscript ^. 



H for the case of two leads, Left (L) and Right (R) con- 
nected to a central region (C), is 





"Hll 


Hlc 







H = 


HcL 


Hcc 


HcR 


(10) 
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where the different block matrices describe the projec- 
tions of the one-body part H of the Hamiltonian onto 
different subregions. They are explicitly given as 

(H„„),. (z) = + %(C/a(2) - m)] 6..', (11) 



(Hcc),^ (z) = [Kj{z) ~ S.jfl] Sa 



(12) 
(13) 



We focus on the dynamical processes occuring in the cen- 
tral region. These are described by the Green's func- 
tion ^cc projected onto region C. We therefore want to 
extract from the block matrix structure for the Green's 
function 



Qll Qlc Glr 

QcL 0CC OcR 
QrL GrC 0RR 



(14) 



an equation for ^cc- The many-body self-energy in Eq. 
([7]) has nonvanishing entries only for indices in region C. 
This is an immediate consequence of the fact that the di- 
agrammatic expansion of the self-energy starts and ends 
with and interaction line which in our case is confined in 
the central region (see last term of Eq. This also 

implies that S^^[5cc] is a functional of Gcc only. From 
these considerations it follows that in the one-particle ba- 
sis the matrix structure of is given as 



->MB 





^cc [yccj 




(15) 



The projection of the equation of motion ([7]) onto regions 
CC and aC yields 

{0,1 - Hcc(^)}ecc(^, z') ^ S{z, z')l + 

r (16) 

5^Hcaeac(^,^')+ / dz-E^^iz,z)gcciz,z') 

for the central region and 

- H„„(z)}e„c(^, z') = Uo,c0cg{z, z') (17) 

for the projection on aC. The latter equation can be 
solved for QaC-, taking into account the boundary condi- 
tions of Eq. dni) , to yield 



Gaciz,z') = J dzgaa{z,z) HacG Cc{z , z') , 



(18) 



where the integral is along the Keldysh contour. Here we 
defined g^^ as the solution of 



jia^l - Hcja(z)|g„„(z, 2') = S{z, z')l 



(19) 



with boundary conditions Eq. The function g^^ 

is the Green's function of the isolated and biased a- 
lead. We wish to stress that a Green's function g^^ 
with boundary conditions Eq. ^ automatically ensures 
the correct boundary conditions for the Qaciz, z') in Eq. 
psp . Any other boundary conditions would not only lead 
to an unphysical transient behavior but also to differ- 
ent steady state results<^ This is the case for, e.g., ini- 
tially uncontacted Hamiltonians in which the equilibrium 
chemical potential of the leads is replaced by the electro- 
chemical potential, i.e., the sum of the chemical potential 
and the bias. 

Taking into account Eq. the first term on the 

righthand side of Eq. becomes 

^Hcaeac(z,z') = f dz-E,^iz,z)gcciz,z'), (20) 
where we have introduced the embedding self-energy 

'Esmiz, z') = ^ Scm,a(-Z, z') = ^ Hca gaa{z, z')'S-aC, 

(21) 

which accounts for the tunneling of electrons from the 
central region to the leads and vice versa. The embed- 
ding self-energies Som.a are independent of the electronic 
interactions and hence of Qoc-, and are therefore com- 
pletely known once the lead Hamiltonians Ha of Eq. ([3]) 
arc specified. Inserting (PO)) back to then gives the 
equation of motion 



{z9,l-Hcc(^)}ecc(^,^') 



(5(2, z')\ + dz 



^CC 



{z, z) Qcq[z, z'). 



(22) 



An adjoint equation can similarly be derived from Eq. 
([5]). Equation is an exact equation for the Green's 
function Qcc- for the class of Hamiltonians of Eq. ([1]), 
provided that an exact expression for [^cc] ^.s a 
functional of Qcc is inserted. In practical implemen- 
tations Eq. ([22|) is converted to a set of coupled real- 
time equations, known as the Kadanoff-Baym equations 
(see Appendix [X)) . These equations are solved by means 
of time-propagation techniques!^ For the case of un- 
perturbed systems the contributions of the integral in 
Eq. (|22p coming from the real-time branches of the 
contour cancel and the integral needs only to be taken 
on the imaginary vertical track. The equation for the 
Green's function then becomes equivalent to the one of 
the equilibrium finite-temperature formalism. In a time- 
dependent situation the vertical track therefore accounts 
for initial correlations due to both many-body interac- 
tions, incorporated in S^c, and contacts with the leads, 
incorporated in Som- In our implementation (see Ap- 
pendix [XJ we always solve the contacted and correlated 
equation first on the the imaginary track, before we prop- 
agate the Green's function in time in the presence of an 
external field. However, to study initial correlations we 
are free to set the embedding and many-body self-energy 
to zero before time-propagation, which is equivalent to 
neglect the vertical track of the contouri^ This would 
correspond to starting with an equilibrium configuration 
that describes an initially uncontacted and noninteract- 
ing central region. This class of initial configurations is 
commonly used in quantum transport calculations, where 
both the interactions and the couplings are considered to 
be switched on in the distant past. The assumption is 
then made that the system thermalizcs before the bias is 
switched on. Even when this assumption is fulfilled there 
are practical difficulties to study transient phenomena, as 
one has to propagate the system until it has thermalized 
before a bias can be switched on. It is therefore an ad- 
vantage of our approach that thermalization assumptions 
are not necessary. 

To solve the equation of motion Eq. we need to 

find an approximation for the many-body self-energy 
S'^^[5cc] as a functional of the Green's function Qcc- 
This approximation can be constructed using diagram- 
matic techniques based on Wick's theorem familiar from 
equilibrium theory^! which can be straightforwardly 
be extended to the case of contour-ordered Green's 
functions.— In our case the perturbative expansion is in 
powers of the two-body interaction and the unperturbed 
system consists of the noninteracting, but contacted and 
biased system. We stress, however, that eventually all 
our expressions are given in terms of fully dressed Green's 
functions leading to fully self-consistent equations for the 
Green's function. This full self-consistency is essential to 
guarantee the satisfaction of the charge conservation law, 
as is discussed in the next section. 



9. 




FIG. 3: Diagrammatic representation of the many-body ap- 
proximations for Sqq . 



C. Charge conservation 

The approximations for S^(?[^cc] that we use in this 
work involve the Hartree-Fock, second Born and GW 
approximation, which are discussed in detail in Refs. 
[ssHsgIIgsIIbsI and are displayed pictorially in Fig. [S] These 
are all examples of so-called conserving approximations 
for the self-energy, that guarantee satisfaction of funda- 
mental conservation laws such as charge conservation. As 
shown by Baym^^, a self-energy approximation is con- 
serving whenever it can be written as the derivative of a 
hmctional <&, i.e. 



-,MB 
■'CCrs 



[Qcc]{z,z') = 



S<S>[Q 



ccj 



SQ. 



CC.sr 



(23) 



This form of the self-energy is by itself not sufficient to 
guarantee that the conservation laws are obeyed. A sec- 
ond condition is that the equations of motion for the 
Green's function need to be solved fully self-consistently 
for this form of the self-energy (see, e.g., Ref. Issl ). For an 
open system, like our central region, charge conservation 
does not imply that the time derivative of the number of 
particles Nc(t) is constant in time. It rather implies that 
the time-derivative of Nc (t) , also known as the displace- 
ment current^ is equal to the sum of the currents that 
flow into the leads. Below we give a proof in which the 
importance of the 4>-derivability is clarified. We start by 
writing the number of particles Nc{t) as (see Eq. ©) 



Ncit) = -iTi-c [Qccit-,t^ 



(24) 



where the trace is taken over all one-particle indices in 
the central region. Subtracting the equation of motion 



from its adjoint and setting z = t- 
yields 



= t+ then 



dNcjt) 
dt 



-2ReTr, 



c 



where Sec = ^ 



MB 
-'CC 



dz'Eccit-,z)Qcciz,t+) 

'(25) 

n. By similar reasonings we 



can calculate the current flowing across the interface 
between lead a and the central region. The total number 
of particles in lead a is Na = —iTta [Qaa{t-, i+)], where 
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the trace is taken over all one-particle indices in lead a. 
Projecting the equation of motion ([7]) on region aa yields 

idr,Qaa{z,z) = b{z , z')\ ■^Ylaa{z)Q aa{.Z , z) 

+ H„c(2)aca(^,^')- (26) 
Subtracting this equation from its adjoint one finds 



Ia{t) 



dt 



2ReTr, [li^cQca{t^.t+)] 



2ReTrc [gca(t-, i+)H„c] ■ 



(27) 



Substituting in this expression the explicit solution (|18p 
for QaC a-s well as the solution for its adjoint Gca we 
can write the current in terms of the embedding self- 
energy So,n,Q as 



lait) = 2ReTrc 



cm, a 



(28) 



Exploiting this result Eq. ((25|) takes the form 
dNcit) 



dt 



= II + Ir- / dz Tr 



(29) 

Charge conservation implies that the integral in Eq. ([29]) 
vanishes. This is a direct consequence of the invariance of 
the functional $ under gauge transformations. Indeed, 
changing the external potential by an arbitrary purely 
time-dependent function Kr{z) (with the boundary con- 
dition Ar(to) = Ar(to — */3)) changes the Green's function 
according to^ 

gcc,r.[A](^, z) = e^^'-(^)gcc,..(^, z')e~*^='"'\ (30) 

as can be checked directly from the equations of motion 
for the Green's function. From its definition Eq. ((23)) it 
follows that the ^-functional consists of closed diagrams 
in terms of the Green's function Qcc- The phase fac- 
tors of Eq. ([50]) thus cancel each other at every vertex 
and therefore $ is independent of the functions A^. This 
implies that 







<5A,(z) 



( dzdz' 



qrs^C 



<5$ 



SQcC.sriz', Z) SAq{z) 



= J: /dz-rfz-'E-B,,(.-.-')^%^#^,(31) 



qrs^C 



5Kq{z) 



where the sums run over all one-particle indices in the 
central region. Here we explicitly used the <I>-derivability 
condition of the self-energy of Eq. . If we now insert 
the derivative of the Green's function with respect to A,, 
from Eq. ([50)) in Eq. ([51]) and evaluate the resulting 



expression in z = t± we obtain the integral in Eq. ([29p . 
Therefore the last term in Eq. ([^5]) vanishes and the time- 
derivative of the number of particles A^c (t) in the central 
region is equal to the sum of the currents that flow into 
the leads. We mention that in the long time limit the 
number of particles in region C is constant provided that 
the system attains a steady state. In this case /l-I-/r ~ 
and we recover the result of Ref. as a special case. 



D. Equation for the time-dependent current 



The time-dependent current in Eq. ()28p accounts for 
the initial many-body and embedding effects. In the ab- 
sence of an external perturbation lait) = at any time. 
The exact vanishing of the current is guaranteed by the 
contribution of the vertical track in the integral. Discard- 
ing this contribution is equivalent to starting with an 
initially uncorrelated and uncontacted system in which 
case there will be some thermalization time during which 
charge fluctuations will give rise to nonzero transient cur- 
rents. 

Equation ([28]) involves an integral over the Keldysh con- 
tour. Using the extended Langreth theorem^ i^^i^^ for the 
contour of Fig. [5] we can express Ia{t) in terms of real 
time and imaginary time integrals 



Ic{t) = 2ReTrc 



/ die^c(^>i)5:l_„(t,t) 

Jto 

+ [ dtg^cit.t)'^em,ait,t) 
/ dTgl^{t,T)T.[^,^{T,t) 

Jq 



(32) 



where we refer to Appendix |^ for the definition of the 
various superscripts. Equation ()32p provides a general- 
ization of the Meir-Wingreen formula^ to the transient 
time-domain. As anticipated the last term in Eq. (|32|) ex- 
plicitly accounts for the effects of initial correlations and 
initial-state dependence. If one assumes that both depen- 
dences are washed out in the long-time limit (t —^ oo), 
then the last term in Eq. p2p vanishes and we can safely 
take the limit to —^ — oo. Furthermore, if in this limit 
the Green's function becomes a a function of the relative 
times only, i.e., Qcc{i^t') Qcc{i — t'), we can Fourier 
transform with respect to the relative time to obtain the 
Green's function Qcci'^) and the self-energy Hcmi^) in 
frequency or energy space. This is typically the case for 
DC bias voltages where limt^oo Ua{t) = Ua- In terms of 
the Fourier transformed quantities Eq. ([5^ reduces to 
the Meir-Wingreen formula2i for the steady state current 

/oo J 

(33) 



where 



r„(u;) = -2Im{S« (.;)} 



(34) 
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AH = -^[egcH-eccM], (35) 

and where is the Fermi distribution for lead a with 
electrochemical potential fi + Ua. This expression has 
been used recently to perform steady state transport cal- 
culations at GW level i^i^i^ The present formalism al- 
lows for an extension of this work to the time-dependent 
regime. 



E. Electron density in the leads 

In our investigations we arc not only interested in cal- 
culating the density in the central region, but are also 
interested in studying the densities in the leads. In the 
following we will therefore derive an equation from which 
these lead densities can be calculated. If we on the right- 
hand side of Eq. ([26)) insert the adjoint of Eq. (flS]) we 
obtain the expression 

idzQaa{z,z') = 5{z,z')l + llaa{z)Qaaiz,z') 

+ J dz'Ein^a{z,z)g^^{z,z'), (36) 

where we defined the inbedding self-energy as 

Sin,a(2, z') = Hqc5cc(z, zOHcq- (37) 

If we solve Eq. ([55)1 in terms of g^^ and take the time 
arguments at t± wc obtain 

+ j dzdzg^a{t_,z)T,i^^a{z, z)ga^{z, t+). 

(38) 

We see from Eq. ^ that with this equation we can ob- 
tain the spin occupation of orbital i in lead a by taking 
r = s = iaa. The integral in Eq. p8|) is taken along 
the Keldysh contour. In practice we solve the Kadanoff- 
Baym equations for Qcc first. After this we construct 
the inbedding self-energy Sin and calculate the lead den- 
sity from Eq. (|38p converted into real time, using the 
conversion table of Rcf. [s^- 



III. NUMERICAL RESULTS 

In this Section we specialize to central regions con- 
sisting of quantum chains modelled using a tight-binding 
parametrization. We studied the case for which the chain 
extends from site 1 to site 4 and is coupled to a left and 
right two-dimensional reservoirs with 9 transverse chan- 
nels in the left and right leads, as illustrated in Fig. [TJ 
The parameters for the system are chosen as follows. The 
longitudinal and transverse nearest neighbor hoppings in 
the leads arc set to = = —2.0, a = L,K, whereas 



the on-site energy a" is set equal to the chemical poten- 
tial, i.e., a" = /i. The leads are therefore half-filled. Pre- 
cise definitions of these parameters can be found in Ap- 
pendix |B1 The endsites of the central chain are coupled 
only to the terminal sites of the central row in both leads 
and the hopping parameters arc V1.5L ~ V4.5R = —0.5 
(sec Appendix |B] for the labeling) . The central chain has 
on-site energies ha = and hoppings hij = — 1.0 between 
neighboring sites i and j. The electron-electron interac- 
tion in the central region has the form Vijki = Vij SuSjk 
with 

^y = i ... „• (39) 

and interaction strength vu = 1.5. For these parameters 
the equilibrium Hartree-Fock levels of the isolated chain 
lie at ei = 0.39, €2 = 1.32, £3 = 3.19, €4 = 4.46. In all our 
simulations the chemical potential is fixed between the 
highest occupied molecular orbital (HOMO) £2 and the 
lowest unoccupied molecular orbital (LUMO) £3 levels at 
fi = 2.26 and the inverse temperature (3 is set to /? = 
90 which corresponds to the zero temperature limit (i.e. 
results do not change anymore for higher values of /?). In 
this work we will consider the case of a suddenly applied 
constant bias at an initial time to, i.e. we take Ua{t) = 
Ua for t > to and Ua{t) = for t < io- Additionally, the 
bias voltage is applied symmetrically to the leads, i.e., 
Ul = —Ur = U, and the total potential drop is 2U. 



A. Keldysh Green's functions in the double-time 
plane 

All physical quantities calculated in our work have 
been extracted from the different components of the 
Keldysh Green's function. Due to their importance we 
decided to present the behavior of the lesser Green's func- 
tion as well as of the right Green's function Q'^ in the 
double-time plane for the Hartrec-Fock approximation. 
The Green's functions corresponding to the 2B and GW 
arc qualitatively similar but show more strongly damped 
oscillations. In Fig. |4] wc display the imaginary part 
of 

^CCMH^^T^') basis of the initial Hartree-Fock 

molecular orbitals, for an applied bias U = 1.2. This ma- 
trix element corresponds to the HOMO level of the molec- 
ular chain. The value of the Green's function on the time 
diagonal, i.e., n/f(i) = Im[^cc jj^(t, t)] gives the level 
occupation number per spin. We see that nnit) decays 
from a value of 1.0 at the initial time to a value of 0.5 at 
time t = 30. An analysis of the LUMO level occupation 
nL{t) shows that almost all the charge is transferred to 
this level. The discharging of the HOMO level and the 
charging of the LUMO level is also clearly observable in 
the dipole moment as it causes a density oscillation in the 
system (see Section llll CI) . When we move away from the 
time-diagonal we consider the time-propagation of holes 
in the HOMO level. We observe a damped oscillation the 
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frequency of which corresponds to the removal energy of 
an electron from the HOMO level, leading to a distinct 
peak in the spectral function (see Section IIII Bl below) . 

The imaginary part of jjj^it, r) within the HF ap- 
proximation is displayed in Fig. [5] for real times between 
t — and t = 30 and imaginary times from t = to 
T = 5. This mixed-time Green's function accounts for 
initial correlations as well as initial embedding effects 
(within the HF approximation only the latter) . At t = 
we have the ground-state Matsubara Green's function 
and as the real time t increases all elements of Q^^^it, t) 
approach zero independently of the value of r. This be- 
havior indicates that initial effects die out in the long- 
time limit and that the decay rate is directly related to 
the time for reaching a steady state. A very similar be- 
havior is found within the 2B and GW approximation 
but with a stronger damping of the oscillations. 



B. Time-dependent current 

The time-dependent current at the right interface be- 
tween the chain and the two-dimensional lead is shown 
in Fig. [HI for the HF, 2B and GW approximations for 
two different values of the applied bias U = 0.8 (weak) 
and 1.2 (strong). The first remarkable feature is that 
the 2B and GW results are in excellent agreement at 
all times both in the weak and strong bias regime while 
the HF current deviates from the correlated results al- 
ready after few time units. This result indicates that 
a chain of 4 atoms is already long enough for screening 
effects to play a crucial role. The 2B and GW approxi- 
mations have in common the first three diagrams of the 
pcrturbative expansion of the many-body self-energy il- 
lustrated in Fig. [31 We thus conclude that the first order 
exchange diagram (Fock) with an interaction screened 




30 



FIG. 4: The imaginary part of the lesser Green's function 
Qqq jj£f(ii,i2) of the central region in molecular orbital basis 
corresponding to the HOMO level of the central chain. Bias 
voltage U = 1.2, HF approximation. 




FIG. 5: The imaginary part of the mixed Green's function 
G'qc Hui^' ''') '^^ central region in molecular orbital basis. 
Bias voltage U — 1.2, HF approximation. 




t 



FIG. 6: Transient currents flowing into the right lead for the 
HF, 2B and GW approximations with the applied bias U = 
0.8 (three lowest curves) and U = 1.2. 

by an electron-hole propagator with a single polarization 
bubble (with fully dressed Green's functions) contains 
the essential physics of the problem. We also wish to em- 
phasize that the 2B approximation includes the so called 
second-order exchange diagram which is also quadratic 
in the interaction. This diagram is less relevant due to 
the restricted phase-space that two electrons in the chain 
have to scatter and exchange. 

We then turn our attention to the spectral function 
which is defined as 

A(T,c^) = -ImTrc J ^e-^[g>c-g<^]{T +\,T -*-). 

(40) 

For values of T after the transients have died out the 
spectral function becomes independent of T. For such 
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times we denote the spectral function by A{uj) and it 
is easy to show that A{uj) = Trc[A(a;)] where A(a;) is 
defined in Eq. ([35|) . This function displays peaks that 
correspond to removal energies (below the chemical po- 
tential) and electron addition energies (above the chem- 
ical potential). The spectral functions of our system are 
displayed in Fig. El At weak bias the HOMO-LUMO 
gap in the HF approximation is fairly the same as the 
equilibrium gap whereas the 2B and GW gaps collapse 
causing both the HOMO and the LUMO to move in the 
bias window. As a consequence the steady-state HF cur- 
rent is notably smaller than the 2B and GW currents. 
This effect has been previously observed by Thygesen^ 
and is confirmed by our time-dependent simulations. 

A new scenario docs, however, emerge in the strong 
bias regime. The HF HOMO and LUMO levels move 
into the bias window and lift the steady-state current 
above the corresponding 2B and GW values. This can 
be explained by observing that the peaks of the HF spec- 
tral function A{uj) are very sharp compared to the rather 
broadened structures in the 2B and GW approximations, 
see Fig. [T] In the correlated case the HOMO and LUMO 
levels can be exploited only partially by the electrons to 
scatter from left to right and we thus observe a suppres- 
sion of the current with respect to the HF case. From 
a mathematical point of view the steady-state current is 
roughly proportional to the integral of A(u;) over the bias 
window which is larger in the HF approximation. 

The time-evolution of the spectral function A{T, co) as 
a function of T is illustrated in Fig. [5] for the case of 
the HF and the 2B approximation. For these results, the 
ground state system was propagated without bias up to 
T = 40 after which a bias was suddenly turned on. The 
HF peaks remain rather sharp during the entire evolution 
and the HOMO-LUMO levels come nearer to each other 
at a constant speed. On the contrary, the broadening of 
the 2B peaks remains small during the initial transient 
regime (up to T = 70) to then increase dramatically. This 
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FIG. 8; Real-time evolution of the spectral function A{T,u}) 
for the HF (left panel) and the 2B approximation (right panel) 
for an applied bias oiU — 1.2. On the horizontal axis the time 
T and the vertical a:x;is the frequency cj. 
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FIG. 7: Spectral functions A{uj) for HF (uppermost plot), 2B 
(middle plot) and GW (bottom plot) approximation with the 
applied bias U = 0.8 (solid line) and U — 1.2 (dashed line). 



FIG. 9: Transient right current Iii{U, t) as a function of ap- 
plied bias voltage and time in the HF (left panel) and 2B 
(right panel) approximations. 



behavior indicates that there is a critical charging time 
after which an enhanced renormalization of quasiparticle 
states takes place causing a substantial reshaping of the 
equilibrium spectral function. 

The time-dependent current at the right interface as a 
function of applied voltage and time is shown in Fig. [9] 
for the HF and 2B approximation. The figures nicely il- 
lustrate how steady state results are obtained from time- 
dependent calculations: after the transients have died 
out we see the formation of the characteristic I-V curves 
familiar from steady state transport calculations. In the 
HF approximation one clearly observes the typical stair- 
case structure with steps that correspond to an applied 
voltage that includes one more resonance in the bias win- 
dow. These steps appear at bias voltages U = 0.9 and 
U = 1.8. This result is corroborated by the left panel 
of Fig. [TO] in which we display the bias-dependent spec- 
tral function for the HF approximation. Here we see a 
sudden shift in the spectral peaks at these voltages. The 
HF results thus bear a close resemblance to the standard 
non-interacting results, the main difference being that 
the HF position of the levels gets renormalized by the 
applied bias. 

We now turn our attention to the 2B approximation 
in the right panel of Fig[9l We notice a clear step at 
bias voltage of U ~ 0.7 but the broadening of the level 
peaks due to quasiparticle collisions completely smears 
out the second step and the current increases smoothly 
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Bias 



FIG. 10: Spectral function A{lo) for the HF (left panel) and 
2B (right panel) approximation, as a function of the bias volt- 
age. For the 2B approximation the spectral functions for bias 
voltages until U = 0.6 were divided by a factor 30 (blue lines 
in the figure) 



as a function of the applied voltage. This is again cor- 
roborated in the right panel of FigdO] where we observe 
a sudden broadening of the spectral function at a bias of 
U = 0.7. To make this effect clearly visible in the figure 
we divided the spectral functions for biases up to ?7 = 0.6 
by a factor of 30. We further notice that for the 2B ap- 
proximation there is a faster gap closing as a function of 
the bias voltage as compared to the HF approximation. 
Very similar results are obtained within the GW approx- 
imation. We can therefore conclude that electronic cor- 
relations beyond Hartree-Fock level have a major impact 
on both transient and steady-state currents. 



C. Time-dependent dipole moment 

To study how the charge redistribute along the chain 
after a bias voltage is switched on we calculated the time- 
dependent dipole moment 



(41) 



sitions between the ground state levels of the central re- 
gion and the electrochemical potentials of the left and 
right leads. However, the oscillations are visible in all ob- 
servable quantities through the oscillations of the Green's 
function discussed in Section [ill Al Detailed information 
on the electronic level structure of the chain can be ex- 
tracted from the Fourier transform of d{t), sec inset in 
Fig. 111! One clearly recognize the presence of sharp 
peaks superimposed to a broad continuum. The peaks 
occur at energies corresponding to electronic transitions 
from lead states at the left/right electrochemical poten- 
tial to chain eigenstates or to intrachain transitions. We 
will denote a transition energy between leads L and R 
and chain eigenstate i by Aeli and AeiR. Similarly we 
will denote a transition energy between states in the cen- 
tral region as Aei^. In the inset of Fig. [11] the main peak 
structures are labeled from the highest to the lowest tran- 
sition energies with letters (a) to (e) and we will use these 
labels to denote the various transitions discussed below. 
The possible transition energies can be determined form 
the position of the peaks in the spectral functions and 
the lead levels. As expected the dominant peak occurs 
at the intrachain transition energy Ae23 ~ 1-5 (c). This 
roughly corresponds to the average of the equilibrium and 
nonequilibrium gaps and, therefore, must be traced back 
to charge fluctuations between the HOMO and LUMO. 
The other observable transition energies are AeL2 ~ 2.0 
(b), AeL3 ~ 0.5 (e) and AeL4 ~ 1-0 (d) from the left 
lead and AeiR w 0.65 (e), Ae2R « 0.4 (e), Aem « 2.0 
(b) and Ae4R w 3.4 (a) from the right lead. Some of 
the peaks with transition energies close to each other 
(AeL2 & AesR (b) and Aels & AeiR & Ae2R(e)) are 
merged together and broadened. The broadening is not 
only due to embedding and many-body effects but also 
to the dynamical renormalization of the position of the 
energy levels. Further information can be extracted from 
the peak intensities. The peak of the AeL4 (d) transition 
is very strong due to the sharpness of that particular res- 



where the the coordinates of the sites of the chain 

(with a lattice spacing of one) with origin between sites 
2 and 3. As observed in Section fill Al the chain remains 
fairly charge neutral during the entire time evolution. 
However, a charge rearrangement occurs as can be seen 
from Fig. [TT] At [/ = 1.2 both the HOMO and the 
LUMO are inside the bias window, the lowest level re- 
mains below and the highest level above. Electrons in 
the initially populated HOMO then move to the empty 
LUMO and get only partially reflected back. This gener- 
ates damped oscillations with the HOMO-LUMO gap as 
the main frequency, a non- vanishing steady value for the 
LUMO population and a partially filled HOMO. Due to 
the different (odd/even) approximate spatial symmetry 
of the HOMO/LUMO levels a net dipole moment devel- 
ops. 

As we pointed out in a recent Letter;^ the oscilla- 
tions in the transient current reflect the electronic tran- 
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FIG. 11: Dipole moment of the central region as a function of 
time for bias U — 1.2. The inset shows the Fourier transform 
of the dipole moment. 
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onance, see Fig. [71 and its initial low population. On 
the contrary, the transition Acli from the left lead to 
the highly populated level ei is extremely weak due to 
the Pauli blockade and not visible. Correlation effects 
beyond Hartree-Fock theory causes a fast damping of all 
sofar discussed transitions. Only the transitions AeL4 (d) 
and Ae23 (c) arc visible in the Fourier spectrum of the 
2B and GW approximation. 



D. Time dependent screened interaction W 

In Fig. [T2]wc show the trace of the lesser component of 
the time-dependent screened interaction of the GW ap- 
proximation in the double-time plane. This interaction 
is defined as W = v + vPW where P is the full polar- 
ization bubbleS^ (with dressed Green's functions) of the 
connected and correlated system, and gives information 
on the strength and efficiency of the dynamical screening 
of the repulsive interactions. The good agreement be- 
tween the 2B and GW approximations implies that the 
dominant contribution to the screening comes from the 
first bubble diagram, that is « vP^v. From Fig. [T^] 
we see that the trace of the imaginary part of W^{t, t) is 
about 3. Considering that the trace of the instantaneous 
bare interaction u is 6 we conclude that the screening di- 
agrams reduce the magnitude of the repulsion by a factor 
of 2. Another interesting feature of the screened inter- 
action is that it decays rather fast when the separation 
of the time arguments increases. From Fig. [12] we see 
that after a time i « 7 the retarded interaction is neg- 
ligibly small. It is worth noting that such a time scale 
is much smaller than the typical time scales to reach a 
steady state, see Fig. [6l 



E. Time-dependent Friedel oscillations in the leads 




10 10 



FIG. 12: Imaginary part of the trace of the screened interac- 
tion {ti,t2) in the GW approximation. 



We implemented the method described in Section III El 
and based on the inbedding technique to investigate the 
electron dynamics in the leads. This study is of spe- 
cial importance since it challenges one of the main as- 
sumption in quantum transport calculations, i.e., that 
the leads remain in thermal equilibrium during the en- 
tire evolution. 

In Fig. [TSjwe show the evolution of the density in the 
two-dimensional 9-row wide leads (see Fig. [1]) after the 
sudden switch-on of a bias voltage. We display snapshots 
of the lead densities at times t = 0, 1.7, 3.6 and 10 where 
up to 10 layers deep into the leads (where to improve the 
visibility we interpolated the density between the sites). 
Since the atomic wire is connected to the central site it 
acts as an impurity and we see density oscillations in 
the leads following diamond-like pattern. These present 
Friedel oscillations that propagate along preferred direc- 
tions. 

The preferred directions in the density pattern can be 
understood from linear response theory. Given a square 
lattice with nearest neig hbor hopping T = = the 
retarded density response function in Fourier space reads 



x(q,w) 



dk /(eic) - /(ek-Kq) 
(27r)2 LO - ek + Ck+q + ir] 
dii /(£k)(ek - Ck-Hq) 
(27r)2 {uj + inY - (ek - Ek+q)^ 



(42) 



where ek = 2r(cosfcx + cosfcy) is the energy dispersion 
and the integral is done over the first Brillouin zone and 
/ is the Fermi distribution function. At half filling the 
Fermi energy is zero and the Fermi surface is a square 
with vertices in (0, ±7r) and (±7r, 0). The dominant con- 
tribution to the integral comes from the values of k close 
to such vertices where the density of states has van Hove 
singularities. The response function x(q = aQjCj = 0), 
with Q = (tt, tt) the nesting vector, is discontinuous for 
a = 1. Indeed, for every occupied k there exists an a < 1 
such that ek+q = Ck < and the integrand diverges at 
zero frequency. On the other hand for a > 1 the vector 
k + q corresponds to an unoccupied state with energy 
Ek+q > and due to the presence of the Fermi function 
the integrand of Eq. pS]) is well behaved even for uj = Q. 
The discontinuity at Q = (tt, tt) is analogous to the dis- 
continuity at 2kp in the electron gas and leads to the 
Friedel oscillations with diamond symmetry observed in 
Fig. [131 By adding reciprocal lattice vectors we find that 
there are four equivalent directions for these Friedel os- 
cillations given by the vectors Q = ±(7r, ±7r). Each of 
these vectors gives in real space rise to a density change 
of the form (5n(r) ^ e**^ '". Therefore a single impurity in 
a 2D lattice induces a cross-shaped density pattern. Due 
to the fact that in our case the lattice ends at the central 
chain, we only observe two arms of this cross. 

The results of Fig. [131 also allows for testing the as- 
sumption of thermal equilibrium in the leads. The equi- 
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librium density [Top-left panel] is essentially the same 
as its equilibrium bulk value at 0.5. After the switch- 
ing of the bias a density corrugation with the shape of 
a diamond starts to propagate deep into the lead. The 
largest deviation from the bulk density value occurs at 
the corners of the diamond and is about 2% at the junc- 
tion while it reduces to about 1% after 10 layers. We also 
verified that the discrepancy is about 3 times larger for 
leads with only three transverse channels. We conclude 
that the change in the lead density goes like the inverse 
of the cross section. Our results suggests that for a mean 
field description of 2D leads with 9 transverse channels 
it is enough to include few atomic layers for an accurate 
self-consistent time-dependent calculations of the Hartree 
potential. 

IV. CONCLUSIONS 

We proposed a timc-dcpcndcnt many-body approach 
based on the real-time propagation of the KB equa- 
tions to tackle quantum transport problems of corre- 
lated electrons. We proved the continuity equation 
for any ^-derivable self-energy, a fundamental prop- 




FIG. 13: Snapshots of the density in left lead for HF approx- 
imation after the bias U = 1.2 switch-on. On the horizontal 
axes the transverse dimension of the lead (9 rows wide, with 
the site connected to the chain in the center) and 10 layers 
deep. Upper panel left: Initial density. Upper panel right: 
density at time t = 1.7, Lower left panel: density at time 
t = 3.6, Lower right panel: density at time t = 10. The upper 
colorbar refers to the initial density in the upper left panel. 
The lower colorbar refers to the remaining pictures. 



erty in non-equilibrium conditions, and generalize the 
Meir-Wingreen formula to account for initial correlations 
and initial embedding effects. This requires an exten- 
sion of the Keldysh contour with the thermal segment 
(to, to — if3) and the consideration of mixed-time Green's 
functions having one real and one imaginary time argu- 
ment. The Keldysh Green's function in the device region 
Qcc is typically used to calculate currents and densities 
in the device. In this work we also developed an exact 
inbedding scheme to extract from Qcc the TD density in 
the leads. 

The theoretical framework and the implementation 
scheme were tested for one-dimensional wires connected 
to two-dimensional leads using different approximations 
for the many-body self-energy. We found that already for 
4-sites wires screening effects play a crucial role. The 2B 
and GW approximations are in excellent agreement at 
all times for moderate interaction strength (of the same 
order of magnitude of the hopping integrals) while the 
HF approximation tends to deviate from the GW and 
2B results after very short times. These differences were 
related to the sharp peaks of the HF spectral function as 
compared to the rather broad structures observed in 2B 
and GW. Our numerical results indicate that the largest 
part of the correlation effects are well described by the 
first bubble diagram of the self-energy, common to both 
the 2B and GW approximation. The screened interac- 
tion was explicitely calculated in the GW approxima- 
tion showing that the screening reduces the interaction 
strenght by a factor of 2 and that retardation effects are 
absent after a time-scale much shorter than the typical 
transient time-scale. The electron dynamics obtained us- 
ing a correlated self-energy differ from the HF dynam- 
ics in many respects: 1) At moderate bias the HOMO- 
LUMO gap closes while in the HF approximation it re- 
mains fairly constant; 2) The HOMO and LUMO reso- 
nances are rather sharp during the transient time to then 
suddenly broaden when approaching the steady state. 
This indicates the occurrence of an enhanced renormal- 
ization of quasiparticle states. The HF widths instead 
remain unaltered. 3) The transient time in the corre- 
lated case is much shorter than in HF, see Fig. [TTl 

The transient behavior of time-dependent quantities 
like the current and dipole moment exhibit oscillations 
of characteristic frequencies that reflect the underlying 
level-structure of the system. Calculating the ultrafast 
response of the device to an external driving field thus 
constitutes an alternative method to gain insight into 
the quasi-particle positions and life-times out of ec^uilib- 
rium. We performed a discrete Fourier analysis of the TD 
dipole-momcnt in the transient regime and related the 
characteristic frequencies to transitions either between 
different levels of the wire or between the levels of the 
wire and the electrochemical potential of the leads. The 
hight of the peaks in the Fourier transform can be inter- 
preted as the amount of density which oscillate between 
the levels of a given transition. In all approximations 
we found that the density mainly sloshes between the 
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HOMO and the LUMO. 

One of the main assumption in quantum transport cal- 
culations is that the leads remain in thermal equilibrium 
and therefore that the bulk density is not affected by the 
presence of the junction. To investigate this assumption 
we considered two-dimensional leads thus going beyond 
the so called wide-band-limit approximation. By virtue 
of an exact inbedding technique we calculated the lead 
density both in and out equilibrium. In the proximity of 
the junction the density exhibits Friedel-like oscillations 
whose period depend on the value of the Fermi momen- 
tum along the given direction. 

In conclusion the real-time-propagation of the KB 
equations for open and inhomogeneous systems provide 
a very powerful tool to study the electron dynamics of a 
typical quantum transport set-up. In this work we con- 
sidered only DC biases. However, more complicated driv- 
ing fields like AC biases or pumping fields can be dealt 
with at the same computational cost and the results will 
be the subject of a future publication. Besides currents 
and densities the MBPT framework also allows for cal- 
culating higher order correlators. It is our intention to 
use the KB equations to study shot-noise in quantum 
junctions using different levels of approximation for the 
Green's function. 



APPENDIX A: THE EMBEDDED 
KADANOFF-BAYM EQUATIONS 



To apply Eq. (|22p in practice we need to transform it 
to real-time equations that we solve by time-propagation. 
This can be done in Eq. (|^^ by considering time- 
arguments of the Green's function and self-energy on dif- 
ferent branches of the contour. We therefore have to 
define these components first. Let us therefore consider 
a function on the Keldysh contour of the general form 

Fiz,z') = F'{z)5(z,z') 

+ e{z, z')F> (z, z') + e{z', z)F<{z, z')ikl) 

where 0{z, z') is a contour Heaviside function^ i.e. 
9{z, z') = 1 for z later than z' on the contour and zero 
otherwise, and 5{z,z') = dz6{z,z') is the contour delta 
function. By restricting the variables z and z' on dif- 
ferent branches of the contour we can define the various 
components of F as 



F^{t,t') 

F\t,T) 

F^{r,t) 



F^\t 



Fit±,to~iT), (A3) 
F{to-iT,t±), (A4) 
-iF{to~iT,to-iT'), (A5) 



and 



F^/^{t,t') = F\t)5{t~t')TB{±tTt')[F> {t,t')- 



(A6) 



For the Green's function there is no singular contribution, 
i.e., = 0, but the self-energy has a singular contribu- 
tion of Hartree-Fock form, i.e., = YF'[g]^ With 
these definitions we can now convert Eq. (|22p to equa- 
tions for the separate components. This is conveniently 
done using the conversion table in Ref. [63- We then 
obtain the following set of equations 



{t,t')+ [sUgr] (^t,t'), 

(A7) 



(t,t') 

(A8) 



idtQ\t,T) = ncc{t)Q\t,T) 



gl - -E' 



+ i 



M 



(r-r'), 



it,r) 

(A9) 
{r,t) 
(AlO) 

(All) 



which are commonly known as the Kadanoff-Baym equa- 
tions. The symbols • and * are a shorthand notation for 
the real-time and imaginary-time convolutions 



[a ■ b] {t,t') 



a{tj)b{t,t')dt, 



(A12) 



[a*b]{t,t') = -i a{t,T)b{T,t')dT. 



In practice we first solve Eq. (jAll[) which describes the 
initial equilibrium Green's function. This equation is de- 
coupled from the other two, since S*^ depends on g^ 
only. The initial conditions for the other Green's func- 
tions and are then determined by g^^ as follows 

a>(o,o) = ig^'{o+), (A13) 



a<(o,o) = ig^\o-), (A14) 

G\0,t) = ig^'{-T), (A15) 

gr(r,0) = tg^'ir). (A16) 

With these initial conditions the Eqs. (|A7[ )- (|A10p can be 
solved using a time-stepping algorithm!^ 

APPENDIX B: EMBEDDING SELF-ENERGY 



From Eq. ((2T|) and Eq. p3)) we see that the embedding 
self-energy has the form 



(Z, Z) = ^ Vfc,ja£/„„^y (Z, z') V,a„ 



(Bl) 
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where k and I label orbitals in the central region. As 
can be seen from this equation, the calculation of the 
embedding self-energy requires the determination ofg^^. 
Since for the isolated lead a the time-dependent field is 
simply a gauge, g^^ is of the form 



9aa{z, z') = g°aa{z, z') exp (^-i J^^ dz Ua{z)^ , 



(B2) 



where g^^ is the Green's function for the unbiased lead, 
and the integral in the exponent is a contour integral. 
The Green's function has the form 

gUz,z')^e{z,z')g''^>{z,z') + 0{z',z)g'^^<{z,z'). (B3) 

It therefore remains to specify gaa- In the following 
we will for convenience separate out the spin part from 



the Green's function and write q5'„ 

act, 





aa.ij ■ 



We will now give give an explicit expression for g^^ 
for the case of two-dimensional leads. The case of three 
dimensions can be treated similarly. We consider a lead 
Hamiltonian of a tight-binding form, that is separable 
in the longitudinal {x) and the transverse (y) directions. 
Therefore the indices in the one-particle matrix hfj of 
Eq. ([3]) denote sites i = {x,y),j = {x',y') where x and 
y are integers running from zero to N" and Ny. At the 
end of the derivation we take the limit N" — > oo. The 
Hamiltonian matrix for the leads is then of the form 



where A and r are matrices that represent longitudinal 
and transverse chains and a" is an on-site energy. Hence 



p 



rat 



(B5) 



where p — {j)x,Py) is a two-dimensional index spanning 
the same one-particle space. The matrix = D'^" (g) 
Z?^" is a direct product of the unitary matrices D'^°' and 
Z?^" that diagonalizc the matrices r" and A" in Eq. (|B4p 
The functions gaa.p have the explicit form 

g°^<Jz,z') = z/(6,„)e-'/^''^-(^--''), (B6) 
9°'>p(^,^') = z(/(ep.)-l)e-^^'^^"(^--^), (B7) 

with /(e) = l/(e^^'~''^ -I- 1) the Fermi distribution func- 
tion. In these expressions epa = ^p^a ~^ ^p^^a' "^here e^^^, 
and e^^^ are the eigenvalues of matrices t" and A". In 
the case the matrices t" and A" represent tight-binding 
chains with nearest neigbour hoppings and and 
zero on-site energy, we have 



D 



X 



2 . nxpoo , 
sm , 

iv^ + 1 ^iv^ + r 

irpx 



2T^ cos( 



NS 



1 



(B8) 
(B9) 




FIG. 14: Tight-binding system for finite 2D leads connected 
to scattering central region. 



and similarly for the transverse transformation matrix 
DJ" and energy el If we insert these expressions in 

yPy ^ Py^ 

Eq. (jB5P and take the limit Nx oo such that we can 
replace summation over px by an integration over the 
angular variable cj) = 7Tpx/{N" 4-1), then wc obtain 



Q^^iAz,z) = > smi ^^isml —) 

y Py—^ ^ ^ 

1 



X — / (i(/)sin(a;0) sin(x'0) 



where now 



2T^ cos( 



2T^ cos ( 



(BIO) 



(Bll) 



(B4) 

The expression for g^^aii obtained from Eq. (jBlOp 



by simply replacing the Fermi function / by / — 1. Let 
us now turn to the embedding self-energy. In this work 
we consider the case that we only have hopping elements 
Vi^ka between central sites i and the first tranverse layer 
of the leads, which are labeled by elements k = {l,y) 
where y = 1 . . . Ny . However, the entire formalism can 
extended to more general cases. This means that wc take 



i.ka 



y%,ya if fc = (1,?/) 







otherwise 



(B12) 



In that case in Eq. (jBlOp only the contribution with 
X = x' = \ survives. Then the product of the sine func- 
tions can be written in terms of the eigenenergies of the 
isolated leads as 



E 

y,y',Py^l 



X sin(^^)sin(^^) 

^ l\Tn/ I 1 ^ ATrv I 1 



V 




X /(e)e-'^-'"*^( 



(B13) 



where we defined En 



e'L ^ and ^ 



e^^^ =b 2|T^|. The expression for S^^.q obtained 
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from Eq. (|B13p by simply replacing the Fermi function / 
by / — 1. In the case that there is no transverse coupling, 
i.e., = 0, the integral is independent of Py and the 



sum over py can be performed to yield Syy/ . Then the 2D 
self-energy becomes a sum of self-energies over separate 
ID leads. 
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